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OH-  THE  STABILITY  OF  FLIGHT  VEHICLES  IN  THE 
LOW  REYNOLDS  NUMBER  NON-LINEAR  REGIME 
by 

Otto  W.  K.  Lee 
Carol  M.  Vaczy 
Eugene  E.  Covert 

Introduction 

v?This  report  addresses  some  Issues  associated  with  dynamic  stability 
of  vehicles  flying  at  speeds  such  that  the  flight  Reynolds  number  is  of 
the  order  of  a  few  hundreds  of  thousands.  In  particular  the  wing 
loading  and  speed  are  such  that  the  flight  occurs  at  relatively  high 
angles  of  attack.  Under  these  circumstances  the  aerodynamic 
characteristics  of  the  vehicle  are  likely  to  be  non-linear  in  nature. 
Here  the  nature  of  the  flight  itself  could  be  significantly  different 
from  that  expected  under  those  conditions  where  the  aerodynamic 
characteristics  are  linear.  Under  those  conditions  the  ideas  underlying 
stability  are  fairly  well  understood.^The  issues  are  discussed  clearly 
in  several  books  that  are  readily  available  to  anyone  interested  in 
learning  more  about  the  topic.  Two  that  come  to  mind  at  once  are 
Etkin's  Airplane  Dynamic  Stability,  and  Control 1  and  McCormack's, 
Aerodynamics,  Aeronautics  and  Flight  Mechanics^.  There  are  other 
excellent  books  and  their  omission  is  not  to  be  interpreted  as  a 
denigration.  It  is  simply  that  the  authors  are  most  familiar  with  those 
two.  Both  these  books  have  a  fairly  complete  history  of  the  study  of 
dynamic  stabiity  of  flying  machines,  starting  with  E.B.  Wilson  in  the 
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United  States,  and  Bryan  in  Great  Britain.  The  subject  was  well 
understood  by  the  late  20' s  and  early  30* s  as  the  reader  may  decide  by 
reading  B.M.  Jones's  article  in  Volume  V  of  Durand's  Aerodynamic 
Theory^.  The  basic  problems  of  the  coupling  of  the  aerodynamic 
characteristics  to  the  motion  is  discussed  there  for  a  variety  of 
special  cases  that  lead  to  understanding  of  the  complete  problem.  The 
modern  reader  may  find  the  presentation  itself  difficult  to  follow 
because  the  current  notation  is  clearer,  and  because  the  examples  are 
somewhat  archaic. 

Indeed  many  improvements  have  been  made  in  the  matter  of  computing 
the  values  of  the  aerodynamic  derivatives  themselves,  as  well  as  in  the 
matter  of  displaying  the  stability  itself.  In  the  present  case  a 
similar  problem  exists,  even  in  comparison  with  the  References  1  and  2. 
There  are  two  reasons  for  this  problem.  First,  in  the  low  Reynolds 
number  regime,  a  number  of  non-linearities  appear  in  the  aerodynamic 
characteristics  which  preclude,  at  least  initially,  the  use  of  the 
so-called  stability  derivatives.  Second,  the  separation  into  the 
symmetric  and  antisymmetric  motion  may  not  be  realistic.  Finally,  there 
is  the  matter  of  determining  the  aerodynamic  characteristics  themselves 
in  the  low  Reynolds  number  regime.  In  the  discussion  that  follows,  each 
of  these  issues  will  be  treated.  The  method  of  approaching  these  issues 
will  follow  the  currently  accepted  spirit  of  mathematical  modeling. 

That  is,  the  goal  is  to  retain  the  primary  physical  processes  in  a 
convenient  form  that  allows  the  results  to  be  representative  of  those 


that  could  be  obtained  if  a  truer  representation  of  the  physical 
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phenomena  either  were  available.  An  example  of  the  latter  might  be  to 
couple  the  rigid  body  motion  to  the  Navier-Stokes  equations,  and  to  a 
representation  of  the  vehicle's  configuration.  In  principle  this  would 
allow  for  a  complete  solution  of  the  motion  along  a  flight  path.  In 
practice,  it  is  not  clear  that  such  a  step  would  be  possible,  let  alone 
practical.  And  even  if  it  were,  the  issue  of  determining  the  stability 
in  a  practical  form  for  discussion  would  still  have  to  be  faced.  In 
this  discussion  an  approximate  way  of  estimating  the  aerodynamic 
characteristic  will  be  discussed  first.  The  issue  of  stability  and  its 
representation  will  be  discussed. 

Equations  of  Motion 

The  equations  of  motion  will  be  written  in  a  standard  form.  The 
origin  of  the  coordinate  system  is  located  at  the  center  of  mass  of  the 
vehicle  and  is  assumed  to  be  fixed.  The  x-axis  points  forward.  The 
y-axis  is  to  the  right  if  one  is  looking  in  the  positive  x  direction. 
The  z-axis  completes  the  orthogonal  set,  and  is  therefore  positive 
downward  from  the  x-y  plane.  The  components  of  linear  velocity  along 
the  x,y,z  axes  are  u,v,w  respectively.  Corresponding  to  the  linear 
velocity  vector,  the  components  of  the  angular  velocity  vector  are 
p,q,r  along  the  x,y,z  axes.  The  components  of  angular  momentum  can  be 
written  as  »  Iij  p^.  in  this  shorthand  the  summation  notation 
is  used  and  p^  corresponds  to  p,  p2  =*  q,  and  p3  =  r.  Note  the 
usual  assumption  of  lateral  symmetry  has  been  made  so  the  products  of 
inertia 

*12  ”  I xy  *  0  and 
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LyZ 


The  product  of  inertia  in  the  pitch  plane 
*13  -  Ixz 

is  normally  small,  and  in  this  case  it  seems  to  be  slightly  negative. 
Thus,  the  yawing  angular  velocity  reduces  the  angular  momentum 
in  roll,  and  conversely.  Ultimately  this  effect  was  neglected  due  to  a 
lack  of  information  about  the  inner  disposition  of  dense  equipment 
inside  the  vehicle.  We  chose  to  write  the  dynamical  equations  in  first 
order  form  as  shown  below: 
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This  set  of  equations  needs  to  be  related  to  a  space  fixed  set  if 
altitude  and  distance  are  important  variables.  This  relation  is 
accomplished  by  defining  a  rotation  matrix  that  puts  the  velocity  vector 
(u,v,w)  into  and  earth  fixed  system  (U,V,W).  Then  one  integrates  the 
resulting  equation  to  give  the  position  vector  (X,Y,Z).  Conceptually, 
and  computationally,  the  easiest  way  to  carry  derive  this  set  is  to 
compute  the  rate  change  of  the  direction  cosines  of  the  X,Y,Z  set  in  the 
x,y,z  frame.  Physically  the  Euler  angles  are  easier  to  visualize,  and 
thus  this  form  was  used,  namely 

$  *»  p  +  q  sec9  tane  sin®  +  r  tans  cos® 
h  =  q  cos®  -  r  sin® 


-q  sin<t>  +  r  sece  cost 


and  the  R  matrix  used  for  the  transformation 


[iO  ■  [E]  •  [i] 


and/ 

~cosi|jcos6  -sinijicos<t>  +  cos\jisin0  sin$  sin\|;siniti  +  cosijjsine  cos$ 

[R]  =*  l  sini|icos6  cos<|jcos$  +  s  initialising  -sin\iisinij>  +  costy sine cos4> | 
_-sine  cos0sin<t>  cosecosti 

Note  here  the  order  of  the  rotation  is  to  yaw  nose  to  the  right  about 
the  body  fixed  z-axis  (or  the  earth  fixed  Z-axis  since  they  are  colinear 
at  this  point)  to  the  angle  ti>.  Then  one  pitches  nose  up  about  the  body 
fixed  y-axis  to  the  angle  0  and  finally  rolls  about  the  body  fixed 
x-axia  to  the  angle  $ .  Note  the  calculation  of  the  Euler  angles 
involves  knowing  the  angular  velocities  about  the  vehicle  body  axis. 

Any  integration  of  these  equations  requires  initial  conditions. 

For  our  purposes  we  have  selected  equilibrium  conditions  followed  by  an 
upset  in  the  velocity  conditions.  We  have  based  the  equilibirum  state 
upon  the  aerodynamic  characteristics  measured  in  the  wind  tunnel7. 
Aerodynamic  Characteristics 

The  calculation  of  the  aerodynamic  characteristics  of  the  vehicle 
may  be  determined  at  several  levels  of  sophistication.  For  the  purposes 
of  this  study  it  is  sufficient  to  assume  that  the  two-dimensional 
airfoil  characteristics,  i.e.,  the  section  lift,  drag,  and  moment 


coefficients,  are  available  as  functions  of  the  angle  of  attack,  at  the 
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flight  Reynolds  number.  It  will  be  assumed  that  thiB  section  data  is 
applicable  to  the  flight  circumstances.  That  is,  it  will  be  assumed  the 
effects  of  vibration  on  the  aerodynamic  coefficients,  either  in  the  wind 
tunnel  or  in  flight,  are  the  same.  It  will  be  further  assumed  that  the 
surface  conditions  are  similar  too.  These  assumptions  are  made 
primarily  in  the  hope  that  the  stability  phenomena  will  be  the  same. 

One  further  assumption  is  necessary.  That  is  the  assumption  of 
quasi-steady  flow.  Undoubtedly  this  assumption  may  not  be  valid,  but 
lacking  data,  or  realistic  computations,  it  is  necessary  to  make  this 
assumption. 

The  simplest  approach  to  the  calculation  of  the  aerodynamic 
characteristics  is  to  use  "strip  theory.”  In  this  approximation  the 
local  sidewash,  upwash,  and  longitudinal  velocity  components  are  used 
to  determine  the  local  flow  angles  and  the  local  dynamic  pressure. 

This  allows  one  to  determine  the  local  section  characteristics  of  the 
wing  and  the  tail.  By  the  appropriate  integration  one  can  determine 
the  instantaneous  forces  and  moments.  In  so  far  as  the  local 
velocities  include  the  effects  of  the  angular  motion  as  well  as  the 
linear  motion,  the  forces  and  moments  reflect  those  effects.  Indeed, 
if  the  rolling  motion  induces  a  stall,  the  change  in  the  pitch  plane 
forces  and  moments  appears  in  a  natural  wa> .  There  are  two  serious 
objections  to  the  use  of  strip  theory.  The  loading  is  not  diminished 
at  the  tip  due  to  the  shed  vorticity,  and  hence  the  tip  effects  are 
weighted  excessively.  This  is  clearly  shown  by  examination  of  the 
results  summarized  by  Betz4.  The  other  serious  objection  is  due  to 
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the  implicit  assumption  that  the  flow  remains  more  or  less  in  the 
relative  free  stream  direction.  Winkleman's  surface  flow  pictures^ 
clearly  show  that  this  implicit  assumption  is  wrong.  However,  until 
better  means  are  available  to  correct  for  this  lack  of  tools,  one  must 
hope  the  net  results  will  not  be  too  greatly  in  error. 

The  limitations  of  strip  theory  can  be  at  least  partially 
alleviated  through  the  use  of  Prandtl’s  integral  equation.  This 
equation  allows  one  to  write  the  distribution  of  downwash  due  to  the 
non-uniform  circulation  distribution  through  an  integral  of  the  spanwise 
drivative  of  the  circulation.  This  derivative  is  closely  related  to  the 
shed  vorticity,  which  induces  the  downwash.  The  local  angle  of  attack 
can  be  computed.  Including  the  downwash,  and  the  section  characteristics 
determined.  The  result  is  a  functional  integral  equation  that  can  be 
solved  by  Iteration.  The  solution  of  this  problem  is  fairly  complicated 
since  the  principal  part  of  the  integral  is  needed.  One  can  expand  the 
integrand  in  a  Fourier  Series,  and  use  Glauert's  integral  formula  to 
obtain  a  complicated  functional  equation.  Fortunately  an  approximation 
due  to  Schrenk5  is  usually  accurate,  and  much  simpler  to  use.  This 
approximation  is  based  upon  a  three-step  process.  First,  one  calculates 
the  strip  theory  loading,  and  integrates  it  to  obtain,  say,  the  lift 
coefficient.  Then,  one  finds  the  quarter  or  half  ellipse  that  has  the 
same  area  under  it  as  the  strip  theory  integral,  which  is  really  an 
algebraic  step.  The  estimate  of  the  actual  span-wise  loading  is  then 
the  average  of  the  strip  theory  and  the  ellipse.  In  the  case  of  a 
rolling  or  yawing  wing  the  elliptical  distribution  is  multiplied  by  a 
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correction  term  of  the  form  (1  +  k  times  2y/b)  to  give  the  estimate  of 
asymmetric  effects  on  the  loading.  The  value  of  "It"  can  be  estimated 
using  the  fundamental  integral  equation,  and  linear  aerodynamics.  The 
analysis  describing  this  use  of  Schrenk’s  approximation  will  be  given 
below. 

Increasing  sophistication  would  require,  effectively,  a  correction 
to  the  camber  of  the  airfoil  due  to  the  effects  of  the  flow  field  due 
to  the  adjacent  sections.  In  view  of  the  fact  that  the  initial 
correction  due  to  the  downwash  is  only  handled  approximately,  it  seems 
unreasonable  to  even  attempt  to  apply  the  higher  order  correction. 
Application  of  Schrenk's  Approximation 

The  calculation  that  was  discussed  above  really  involves  two 
steps.  The  first  step  is  the  determination  of  the  local  geometric 
angle  of  attack.  This  follows  from  the  determination  of  the  local 
linear  and  angular  velocities.  in  the  case  of  the  horizontal  tail  one 
must  account  for  the  downwash  due  to  the  forward  wing,  as  well.  After 
the  local  angle  of  attack  is  known,  one  can  follow  Schrenk’s  procedure 
for  all  the  symmetric  cases.  The  procedure  must  be  modified  in 
application  to  the  asymmetric  motions  due  to  rolling  and  yawing.  In 
this  case  an  extra  term,  which  contributes  zero  net  force,  is 
proportional  to 

( 2y/b )  -^1  -  (2y/b)2 

or  sine  two  theta  where  theta  is  arccosine  (2y/b).  The  constant  of 
proportionality  is  found  such  that  the  area  under  the  strip  theory 
loading  on  one  halfBpan  curve  has  the  same  value  as  that  described  by 


the  terms  above.  In  the  particular  geometry  used  here,  a  result  that 
was  essentially  equal  to  that  obtained  by  the  Schrenk  procedure  could  be 
found  by  simply  using  strip  theory  to  calculate  the  lift,  say,  and  then 
comparing  that  estimate  with  the  measured  lift  from  the  wind  tunnel  in 
ratio  form.  The  resulting  ratio  was  then  applied  to  all  strip  theory 
results.  This  step  reduced  the  amount  of  calculation  considerably.  The 
iosa  in  accuracy  was  such  that  the  damping  terms  due  to  angular  motion 
were  too  small,  so  the  procedure  was  felt  to  be  conservative. 

Analyis  of  Stability 

The  analysis  of  stability  will  follow  three  specific  steps.  First, 
the  stability  derivatives  given  in  Ref.  7  will  be  used  in  a  classical 
stability  analysis  following  the  procedures  outlined  by  Etkinl .  The 
elementary  pitch  plane  analysis  will  be  applied  to  the  pitch  plane  only 
that  is  based  upon  the  local  value  of  the  aerodynamic  loads.  Finally, 
an  analysis  of  the  full  non-linear  equations  will  be  presented.  This 
analysis  will  be  based  upon  the  ideas  of  Liapunov^.  Liapunov 
stability  is  based  upon  the  concept  of  the  phase  plane.  The  phase  plane 
is  a  plane  whose  coordinates  are  the  velocity  and  the  displacement.  The 
dynamic  motion  is  represented  as  a  locus,  or  trajectory,  in  the  phase 
plane.  Each  point  on  the  trajectory  corresponds  to  the  instantaneous 
state,  that  is,  the  instantaneous  velocity  and  position.  An  equilibrium 
position  is  a  point  on  this  plane,  if  the  equilibrium  is  a  steady  state. 
If  the  equilibrium  corresponds  to  a  limit  cycle  then  the  trajectory  in 
the  phase  plane  is  a  closed  curve.  Naturally  if  the  motion  traces  out  a 


line  that  moves  further  and  further  from  the  initial  point,  the  motion 
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is  said  to  be  unstable.  Liapunov  has  summarized  this  succinctly, 
namely:  "If  a  point  on  a  stable  trajectory  in  the  phase  plane  is 

perturbed  slightly  and  if  the  resulting  motion  remains  close  to  the 
original  trajectory  for  all  time,  the  dynamical  system  whose  motion  is 
represented  by  the  trajectory  is  said  to  be  stable."  We  chose  to 
represent  the  phase  plane  as  the  total  length  of  the  perturbed  velocity 
and  its  integral  with  respect  to  time  for  the  ordinate  and  the  abcissa, 
respectively.  The  perturbation  is  about  the  equilibrium  state,  which  is 
denoted  by  the  subscript  0.  Thus  the  perturbed  velocity  vector  may  be 
written  as  (u  -  ug,  v,  w  -  wo,  RtP.  Rn<3  ■  Rbr)-  Here  Rt, 

Rn,  and  Rb  are  the  appropriate  components  of  the  torsion,  the  normal 
and  the  binormal  radii  of  curvature  of  the  vehicle's  trajectory  in 
physical  space.  If  one  is  not  interested  in  comparing  the  stability  of 
several  different  systems,  but  rather  the  stability  of  one  system  under 
different  circumstances  such  as  gross  weight,  then  if  the  equilibrium 
state  corresponds  to  steady  motion,  if  it  is  sufficient  to  represent  the 
velocity  vector  as  (u  -  uo ,  v,  w  -  wg ,  p,  q,  r)  since  the  product  of 
the  several  radii  and  the  several  angular  velocities  approach  zero  as 
the  system  approaches  equilibrium.  The  simpler  approach  is  adequate  for 
the  purpose  of  discussing  stability. 

Specific  Application 

The  vehicle  considered  in  this  investigation  was  the  Long  Duration 
Expendable  Decoy  (LODED)  formulated  by  Locus,  Inc.  for  the  Naval 
Research  Laboratory.  Data  for  the  vehicle  and  its  wind  tunnel  tests  at 
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the  University  of  Maryland  were  obtained  from  the  LODED  Summary 


11 


Report7 . 

Physical  dimensions  used  were  based  on  Figs.  1  and  2  and  the 
supplied  airplane  characteristics  data  sheets  (see  Appendix  A).  The 
configuration  selected  was  with  lower  winglets,  wing  incidence  +10  degrees, 
tail  incidence  +5  degrees,  and  center  of  gravity  at  the  balance  center 
of  the  wind  tunnel  model.  The  airfoil  used  for  the  wing,  tail,  and 
winglets  was  the  Wortman  FX63-137.  In  addition  to  this,  the  Lissaman 
airfoil  was  used  to  investigate  the  effect  of  hysteresis  on  stability 
calculations.  (Section  lift  for  both  cases  shown  in  Figs.  3A  and  3B). 

Three  methods  were  used  to  numerically  model  the  dynamic  response 
of  the  aircraft.  The  first  is  a  longitudinal  plane  method  that  only 
solves  for  the  velocities  u,  w,  the  flight  path  angle  y,  the  fuselage 
reference  line  angle  0,  and  the  angle  of  attack  a.  The  method  used  the 
linear  stability  derivatives  supplied  in  a  printout  that  came  with  the 
data.  The  equations  for  this  method  are  different  than  those  for  the 
other  two  methods,  and  will  be  discussed  further  under  the  section 
entitled  Pitch  Plane.  The  second  method  used  is  also  a  linear  method. 

It  solves  for  the  external  forces  and  moments  using  linear  stability 
derivatives,  and  uses  these  in  the  nonlinear  equations  of  motion.  The 
third  method  used  is  a  fully  nonlinear  method.  It  uses  strip  theory  to 
solve  for  the  external  forces  and  moments,  and  then  corrects  this  result 
for  real  effects  using  a  wind  tunnel  correction  factor.  Using  this 
method  the  response  of  the  vehicle  at  high  angle  of  attack  can  be 
modelled.  Figure  3C  shows  the  lift  coefficient  for  the  wind  tunnel 
test,  along  with  the  slope  used  for  the  stability  derivatives.  This  big 


difference  in  slope  (especially  at  the  a=0  point),  the  linear  and 
nonlinear  methods  cannot  be  reliably  compared  if  the  motion  varies  over 
a  wide  a  range. 

The  basic  equations  of  motion  were  described  above  and  come  from 
Ref.  1.  These  equations  have  been  modified  to  accommodate  the  pitch 
roll  and  the  pitch  yaw  products  of  inertia,  even  though  for  this  vehicle 
they  equal  zero.  The  coordinate  system  used  was  body  axis,  shown  in 
Fig.  4.  The  controls  were  assumed  to  be  fixed,  so  that  the  total  number 
of  independent  equations  is  reduced  to  9,  matching  the  number  of 
unknowns.  The  three  equations  to  determine  the  position  of  the  center 
of  gravity  at  every  time  were  also  included,  although  they  were  never 
needed  due  to  a  fixed  density  assumption. 

The  equilibrium  condition  was  assumed  to  be  8=0.  It  was  desired  to 
analyze  the  behavior  at  a  fixed  Reynolds  number  of  285,000.  This  was 
accomplished  by  varying  the  mass  of  the  aircraft  for  different  initial 
lift  coefficients.  The  mass  moments  and  products  of  inertia  were  varied 
only  by  the  mass  as  the  radius  of  gyration  was  assumed  to  be  constant. 
The  mass  moment  and  products  of  inertia  for  the  propeller  were  assumed 
to  be  zero. 

Method  of  Solution 

The  nine  independent  equations  can  be  put  into  a  matrix  form  so 
that  all  of  the  variables  can  be  solved  for  at  the  same  time.  Since 
this  results  in  one  basic  first  order  nonlinear  ordinary  differential 
equation,  the  system  can  be  solved  using  a  marching  method.  The 
initial  conditions  are  supplied,  and  the  new  values  are  calculated  in 
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terms  of  the  old  values.  This  results  in 
U(n+1)  =*  U{n)  +  At  *  6u 

(where  in  this  context  U  represents  the  matrix).  It  was  found  that  a 
time  step  above  0.1  seconds  caused  the  numerical  method  to  be  unstable, 
and  so  all  the  runs  were  made  with  a  time  increment  of  0.1  seconds. 

Stability  Analysis 

Linear  Stability  Analysis 

Linear  stability  analysis  was  taken  from  Ref.  1  for  both  the 
longitudinal  and  lateral  modes.  Both  were  solved  numerically, 
revealing  the  short  period  mode  along  with  the  phugoid  mode.  The 
periods  for  the  modes  were  calculated,  and  are  discussed  in  the  linear 
stability  derivative  section. 

Nonlinear  Stability  Analysis 

Liapunov  Stability  analysis  was  used  to  evaluate  the  stability 
for  both  the  linear  and  nonlinear  methods.  In  the  phase  plane 
(Jt-vs-x),  if  a  point  on  a  trajectory  is  displaced  a  small  amount, 
the  displaced  point  remains  arbitrarily  close  to  the  original 
trajectory.  The  trajectory  will  look  similar  to  the  spiral  shown  in 
Fig.  5.  The  difficulty  in  uniquely  determining  the  sign  of  each 
contribution  is  overcome  by  assuming  continuity  with  increasing  time. 


Thrust  and  Propeller  Effects 


Assuming  the  propeller  is  constant  speed,  the  variation  in  thrust 
depends  on  the  instantaneous  forward  flight  speed  (Vx). 

Thrust  =■  T0.(i  -  (vx-u0)/Uco) 

where  subscript  o  indicates  equilibrium  value;  u  =  u  velocity 
component,  »  total  velocity. 

Forces  are  developed  by  the  propeller  which  are  normal  to  its  axis 
of  rotation.  This  is  due  to  the  section  lift  characteristics  of  the 
propeller  blades.  Thus  a  "lift"  (negative  z-force)  is  generated  if  the 
propeller  is  at  an  angle  of  attack,  aprop.  Similarly,  a  side  force  is 
generated  if  the  propellers  are  in  sideslip,  Bprop.  Since  the 
propeller  is  at  the  tail  of  the  vehicle,  downwash  must  be  accounted  for; 

“prop  “  .5*a  for  the  aircraft.  Based  on  the  usual  approximations  for 
propellors  and  a  typical  section  lift  curve  slope  is  .05  (1/deg),  the 
forces  are: 

F  -  -.5pV2S  (2.865)  a 
zp  M  prop 

F  =■  -.5pV2S  (2.865)  P 
*p  prop 

And  resulting  pitch  and  yawing  moments: 


N  -  F  *  X 

p  yp  p 

where  Xp  »  location  of  propeller  from  c.g. 
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Linearized  Stability  Analysis  of  Controls-Flxed  Motion 

The  dynamic  response  of  the  LODED  vehicle  was  analyzed  using  the 
supplied  lateral  and  longitudinal  derivatives.  The  vehicle  is 
perturbed  from  its  equilibrium  position  by  a  triangular  impulse  to  one 
of  its  velocity  components.  Its  linear  and  angular  velocities  and 
orientation  are  calculated  versus  time.  In  addition,  the  stability 
determinant  is  expanded  and  the  characteristic  roots  are  obtained. 

Equilibrium  lift  coefficients  chosen  were  Cl  -  .6  and  Cl  -  1 . 

For  the  desired  Reynolds  number,  a  lift  coefficient  beyond  this  range 
would  have  resulted  in  unreasonable  thrust  and  attitude  requirements. 
The  angle  of  attack  reported  is  with  respect  to  the  zero  lift  line, 
such  that  at  a  *  0  degrees.  Cl  =  0. 

The  stability  derivatives  were  assumed  to  be  given  in  wind  axes  so 
an  appropriate  axis  transformation  was  implemented.  The  stability 
derivatives  with  respect  to  the  velocity  u  were  calculated,  found  to  be 
small  and  were  neglected. 

Longitudinal 

Twelve  runs  were  made  with  an  input  disturbance  to  the  w  velocity 
component  or  the  u  velocity  component.  The  change  in  the  velocity  was 
such  that  if  all  the  other  motion  variables  were  held  fixed  over  the 
time  of  the  disturbance,  the  desired  change  in  angle  of  attack  (6a) 
would  have  resulted.  Table  I  lists  the  pertinent  information  for  each 
run.  Four  variables  were  plotted  for  each  run,  two  over  a  short  period 

3  sec.)  and  two  over  a  longer  period  (~  30  sec.).  The  long 
period  plots  were  generated  every  second  and  at  each  pitch  rate  (q) 
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zero-crossing  and  not  at  each  time  increment  interval  (.01  sec.).  Thus, 
these  curves  are  not  completely  smooth  and  for  qualitative  purposes 
only.  Runs  for  a  longer  period  (~  200  sec.)  were  to  examine  absolute 
convergence.  Figures  6  and  7  show  typical  results.  The  calculated 
characteristics  of  the  short  period  (A)  and  phugoid  (B)  modes  are: 

Cl  -.6  6o«0  deg. 

(A)  roots:  (-.03375  +  .068281) 

period:  .648  sec.;  time  to-half-amp:  .144  sec. 

un  =  .0762  (1/sec)  ?  -  .443 

(B)  roots:  (-8.4399E-05  +  .00317i) 

period:  13.951  sec.;  time  to-half-amp:  57.574  sec. 

un  -  .00317  (1/sec)  ;  =  .027 

CL  “  1 ■  0o  ”  0  deg. 

(A)  roots:  (-.02027  +  .05485i) 

period:  .807  sec.  time  to-half-amp:  .24  sec. 

un  =  .0585  (1/sec)  x,  =  .347 

(B)  roots:  (-5.7193e-05  +  .00318i) 

period:  13.88  sec.  time  to-half-amp:  84.9  sec. 

un  =  .00319  (1/sec);  5  =  .018 

The  numerical  results  similarly  indicated  dynamic  stability  for  all 
conditions.  As  expected,  the  magnitude  and  duration  of  the  disturbance 
(6a)  affects  only  the  amplitude  of  motion,  not  its  damping  or  frequency 
(period).  There  is  agreement  in  that  the  short  period  and  phugoid 
period  increase  with  a  higher  lift  coefficient,  the  former  more  than  the 


latter.  However,  both  periods  were  shorter  in  the  numerical  solution 
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than  the  periods  obtained  from  the  characteristic  roots.  From  the  runs 
at  C l».6,  the  short  period  ~  0.47  sec.  and  at  Cl”1  ,  the  short 
period  ~  .57  sec.  At  Cl”. 6,  the  phugoid  ~  10.1  sec.,  at  Cl=1, 
the  phugoid  ~  10.08  sec. 

The  differences  may  be  caused  by  employing  a  flight  condition 
different  from  that  at  which  the  stability  derivatives  were  calculated. 
This  includes  the  Reynolds  number,  density,  mass,  and  geometry  (wing 
incidence,  tail  incidence).  The  effect  of  each  of  these  parameters  on 
the  stability  derivatives  must  be  analyzed.  Inclusion  of  the 
previously  neglected  u-derivatives  and  revised  modeling  of  the  thrust 
variation  should  be  considered  in  any  further  runs. 

La teral 

A  perturbation  to  the  v-velocity  component  was  initiated  to  develop 
a  sideslip,  B.  Two  runs  were  performed,  one  at  each  lift  coefficient. 
The  calculated  characteristics  of  the  dutch  roll  (A)  and  spiral 
divergence  (B)  modes  are: 


CL  = 

.  6 ,  9q 

”  0  degrees 

(A) 

root: 

(-.92835  +  Oi) 

time  to-half-amp: 

.052  sec 

root: 

(-.03531  +  .  2Q37i ) 

period 

:  2.17  sec.  time  to-half-amp:  1 

.375  sec. 

un  »  • 

207  (1/sec)  c 

=  .171 

(B) 

root: 

(.0061  +  Oi) 

time  to-double-amp: 

3.01  sec 

- 

f  »  9o 

=  0  degrees 

(A) 

root: 

(-.5589  +  Oi) 

time  to-half-amp: 

.087  sec. 

root : 


(-.0198  +  . 1 60i ) 
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period:  2.77  sec.  time  to-half-amp:  2.45  sec. 

un  =  .161  (1/sec)  ;  =  .1231 

(B)  root:  (.0059  +  Oi)  time  to-double~amp:  8.26  sec. 

It  was  noted  that  for  both  lift  coefficients  the  spiral  mode  is 
divergent.  That  is,  increases  continually  and  is  not  converging  to 
any  particular  heading.  Due  to  the  large  value  of  damping  in  roll,  the 
rolling  mode  converges  quickly  and  is  not  discernible  in  the  results. 

The  period  of  the  "dutch  roll"  mode  from  the  motion  is  ~  0.75  seconds 
for  Cl=.6,  and  ~  .8  seconds  for  0^=1;  again  shorter  than  the 
period  derived  from  the  characteristic  roots  but  still  following  the 
same  trend:  longer  periods  at  higher  lift  coefficients. 

In  both  runs,  the  longitudinal  mode  became  excited  after  a  lateral 
disturbance.  A  cross-coupling  between  the  two  modes  of  motion  arose 
because  the  use  of  the  NACA  stability  axes  for  computing  derivatives. 
Thus  the  angular  velocities  are  transformed  between  body  and  stability 
axes.  This  coupling  is  absent  in  the  longitudinal  disturbances  only 
because  the  reference  (equilibrium)  sideslip  is  zero.  This  is  not  a 
problem  if  one  uses  the  body  axes  consistent  for  the  calculations  of  the 
stability  derivatives.  In  the  linear  regime  this  is  sufficient  to 
"de-couple"  the  modes. 

Pitch  Plane  Solution 

The  second  method  of  solution  used  to  study  the  dynamic  response  of 
the  aircraft  was  one  which  just  considered  the  motion  in  the  pitching 
plane.  The  standard  equations  are  used  (Ref.  1,  pp.  181-185)  except 


that  the  actual  variation  of  with  angle  of  attack  was  used.  They 
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were  solved  by  a  backwards  differencing  technique  which  yields: 

ie  -  g 

0 (n+2 )  =  0  <n+ 1 )  x  2.  -  0(n)  +  At^  G/I 
The  time  step  used  was  0.05  seconds.  For  time  increments  larger  than 
this,  the  numerical  method  was  unstable.  Table  II  shows  the  cases  that 
were  run.  Using  this  method,  the  aircraft  is  stable  for  all 
disturbances  used. 

The  phugoid  period  of  the  calculated  motion  was  ~  22  sec.  This 
was  significantly  longer  than  the  value  calculated  from  the  stability 
derivatives.  All  variables  demonstrated  extreme  sensitivity  to 
disturbances  to  the  pitch  angle  0 . 

The  Nonlinear  Method 
Strip  Theory 

In  order  to  calculate  the  external  forces  and  moments  acting  on 
the  aircraft,  strip  theory  was  used.  This  involves  calculating  the 
local  angle  of  attack  at  each  strip  along  the  surface  (surface  being 
either  wing,  tail  or  winglet),  and  integrating  the  contribution  over 
the  entire  span.  The  section  characteristics  are  assumed  valid  at  each 
strip.  Wind  tunnel  data  and  section  characteristics  were  curve-fitted 
using  IMSL  subroutines  at  the  Joint  Computer  Facility.  Average  error 
for  each  fit  was  below  6%.  Longitudinal  force  and  moment  coefficients 
(Cl,  Cq,  C^)  were  assumed  symmetric  for  positive  and  negative 
sideslip  while  the  lateral  coefficients  (Cy,  Cj_,  Cfj)  were  assumed 
anti-symmetric. 

The  local  velocity  at  each  point  was  given  by  the  following 


/**  • 
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equations, 

ifv=u+uxr 
vx  =  u  +  qz  -  ry 
Vy  =  v  +  rx  -  pz 
vz  =  w  +  py  -  qx 

which  leads  to  local  angle  of  attack,  of: 

a  -  tan"1  vz/vx 
and  aircraft  angle  of  attack: 
aQ  =  tan-1  w/u 

The  velocity  at  each  point  is  now: 

V  -  v£  +  y2  +  v|. 

For  the  wing,  the  local  angle  of  attack  is: 
a  “  tan-1  vz/vx  +  iw 

Using  the  velocity  at  each  point,  the  thrust  now  uses  the  velocity  Vx 
which  accounts  for  the  distance  away  from  the  center  of  gravity. 

The  angle  of  attack  at  the  tail  included  two  correction  factors. 
The  first  takes  care  of  the  time  lag  from  wing  to  tail,  and  the  second 
takes  care  of  the  downwash  at  the  tail.  This  results  in  the  angle  of 
attack  at  the  tail  being  given  by  the  following  equation: 

a  -  1/2  ( tan-1 ( vz/vy ) )  +  it 
The  angle  of  sideslip  is  given  by: 

B  =  tan-1  (vy/vx) 

Since  the  winglets  are  the  same  airfoils  as  the  wing  and  tail,  the 
above  angle  is  used  to  get  a  side  force  from  the  CL-vs-a  curve.  Note 
that  although  the  winglets  are  at  a  different  Reynolds  number  than  the 


wing  and  tail,  insufficient  data  exists  to  correct  for  this  factor. 


In  order  to  correct  strip  theory  for  real  effects  (interference 
with  fuselage,  induced  velocities,  etc.)  a  wind  tunnel  correction 
factor  was  included.  This  factor,  R,  was  equal  to  the  wind  tunnel 
results  divided  by  the  strip  theory  results  if  the  angular  velocities 
are  zero.  The  equations  for  the  external  moments  and  forces  become: 

EXTERNAL  FORCE  (OR  MOMENT)  REAL  = 

EXTERNAL  FORCE  (OR  MOMENT)  STRIP  THEORY  X 
(EXTERNAL  FORCE  (OR  MOMENT)  WIND  TUNNEL)/ 

(EXTERNAL  FORCE  (OR  MOMENT)  STRIP  THEORY  WITH  NO  ANGULAR  VELOCITY) 

The  equations  were  solved  numerically  using  the  trapezoidal  rule. 
The  contribution  to  the  force  or  moment  was  evaluated  at  two  points, 
averaged,  and  then  applied  at  the  strip  between  them.  The  fuselage 
contribution  is  given  by  the  method  previously  discussed  (Appendix  C). 
The  type  of  disturbance  used  is  the  same  as  that  described  in  the  linear 
section. 

Runs  Using  Strip  Theory 

A  summary  of  runs  is  listed  in  Table  III.  Detailed  plots  and 
all  computer  printouts  are  available  upon  request. 

Discussion  of  Strip  Theory  Results 

Typical  plots  are  given  in  Figure  8.  A  few  general  trends  can  be 
identified.  In  the  short  period  mode,  period  about  .1  sec.,  variations 
in  angle  of  attack  are  a  function  of  both  the  disturbance  magnitude  and 
duration.  The  period  is  only  slightly  dependent  on  the  magnitude,  while 
very  dependent  on  the  duration.  For  longer  duration  of  the  disturbance, 
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the  period  is  more  affected  by  the  magnitude  than  in  the  case  of  the 
shorter  duration.  When  the  initial  angle  of  attack  is  in  the  nonlinear 
range,  the  period  becomes  slightly  longer,  on  the  order  of  5%. 

The  damping  in  the  short  period  mode  is  relatively  unaffected  by 
the  magnitude  of  the  disturbance,  although  it  does  seem  more  damped  for 
longer  duration  of  disturbance.  As  expected,  the  damping  in  the 
nonlinear  range  is  less  than  that  in  the  linear  range.  In  general,  in 
both  linear  and  nonlinear  regions,  the  longer  duration  disturbance  is 
more  heavily  damped. 

The  long  period  mode,  the  phugoid,  is  generally  lightly  damped. 

This  is  in  contrast  to  the  short  period  mode  which  is  generally  very 
heavily  damped.  In  the  nonlinear  region,  this  trend  no  longer  is  seen. 
Instead  the  short  period  mode  is  lightly  damped,  while  the  phugoid  is 
the  more  heavily  damped.  Although  the  amplitude  in  the  nonlinear  range 
is  altered  by  the  non-linear  effect,  the  period  is  relatively 
independent  as  expected.  For  both  cases  the  period  phugoid  is  about  10 
seconds,  about  30%  less  than  the  linear  analysis  prediction,  which  is 
frequently  the  case  (Ref.  1). 

Liapunov  stability  indicates  that  the  aircraft  is  unstable, 
primarily  due  to  the  spiral  instability  mode.  Originally  the  phase 
plane  plot  indicates  a  stable  spiral  (Figs.  9-16).  As  time  goes  on,  the 
trajectory  starts  to  move  out  away  from  the  origin.  Since  we  did  not 
suppress  the  lateral  modes,  they  appeared  following  the  disturbance  of 
the  longitudinal  mode.  All  runs  showed  this  instability,  but  to  varying 


degrees.  This  can  be  explained  by  noise  and  round-off  error,  which  is 
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amplified  in  the  unstable  lateral  mode. 

Figures  9,  10,  11  show  the  effect  of  input  disturbance  upon  the 
Liapunov  stability  for  lift  coefficient  of  2.33.  Generally  speaking  the 
results  are  similar  (although  three  different  scales  were  used  to 
illustrate  the  following  patterns).  Figure  9,  the  small  perturbation, 
behaves  nearly  the  same  way  as  the  motion  resulting  from  linear 
aerodynamics.  The  larger  disturbance  pushes  the  aircraft  near  the  stall 
so  the  amplitude  increases  and  the  damping  is  reduced.  Note  some 
saturation  on  the  plotter  is  evident  earlier  in  Figs.  10  and  11.  The 
drift  into  the  lateral  instability  always  occurs  at  about  the  same 
distance . 

Figures  12,  13,  14  show  the  effect  increasing  the  equilibrium  lift 
coefficient  to  2.77.  Here  the  motion  proceeds  past  stall,  and  the 
motion  is  less  oscillatory,  i.e.,  the  "spring  strength"  is  reduced.  In 
fact  the  oscillation  fails  to  build  up  before  the  lateral  instability 
builds  up. 

Figure  15  shows  a  lateral  stability  run,  with  a  lightly  developed 
Dutchroll  and  which  also  shows  strong  spiral  divergence. 

Figure  16  shows  the  hysteresis  loop  (Fig.  3c)  effects  which  shows 
a  pitch  oscillation  that  is  lightly  damped.  It  suggests  that  without 
care,  the  hysteresis  could  be  coupled  with  dynamics  for  a  new  class  of 
instability. 

The  Lissaman  (Fig.  3c)  airfoil  was  used  due  to  an  "interesting" 
hysteresis  loop  in  stall.  Since  the  aircraft  was  not  tested  with  this 


airfoil,  there  is  not  an  appropriate  wind  tunnel  correction  factor.  It 
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was  felt  that  the  correction  factor  for  the  FX-Wortman  airfoil  would  be 
better  than  no  correction  at  all.  The  only  change  made  was  in  the 
incidence  of  the  wing  and  tail.  They  were  changed  at  each  initial  angle 
of  attack  of  the  aircraft  so  that  they  would  reach  the  bottom  loop  at 
the  disturbance.  This  was  done  so  that  the  airfoil  would  be  running  in 
the  range  of  the  hysteresis  loop.  Note  that  the  value  of  mass  and 
equilibrium  thrust  wre  calculated  for  the  FX-Wortman. 

The  Lissaman  airfoil  was  run  under  initial  conditions  that  were 
very  close  to  the  upper  angle  of  attack  end  of  the  lift  curve.  The 
angle  of  attack  was  then  perturbed  in  the  negative  direction,  causing 
the  reaction  to  drive  it  to  the  lower  hysteresis  loop.  Due  to  the 
nature  of  the  algorithm,  once  there  it  was  never  able  to  reach  the  upper 
curve  again.  Hence  this  result  represents  oscillation  into  a  stall 
regime.  Results  for  the  Lissaman  airfoil  are  somewhat  different  from 
those  for  the  Wortman  airfoil.  The  short  period  mode  with  the  Lissaman 
airfoil  is  longer  than  that  with  the  Wortman,  while  the  phugoid  appears 
to  be  slightly  smaller. 

Due  to  lack  of  time,  only  one  lateral  run  was  made.  The  side  slip 
velocity  appeared  to  damp  out  very  rapidly.  The  short  period  mode 
seemed  longer  than  that  in  the  longitudinal  case,  although  it  was  not 
plotted  for  long  enough  to  positively  confirm  this.  The  long  period 
mode  has  a  period  of  about  5  sec. ,  which  is  half  as  short  as  that  in 
the  longitudinal  case.  The  Liapunov  stability  analysis  revealed  a 
slowly  decreasing  velocity  and  a  rather  rapidly  increasing  distance. 


This,  along  with  the  printout  for  the  heading  angle,  seems  to  indicate 


that  the  aircraft  is  seeking  a  new  heading,  i.e.,  the  machine  possesses 
a  spiral,  instability,  or  predicted  by  the  linear  theory. 

Low  Reynolds  Number  Data 

It  is  perhaps  worth  noting  that  aerodynamic  data  that  is  obtained 
at  low  Reynolds  numbers  may  lack  the  precision  and  repeatability  of 
that  taken  at  higher  Reynolds  numbers.  There  are  several  reasons  for 
this  difficulty.  First  the  dynamic  pressures  are  lower  so  the  size  of 
the  forces  and  moments  that  are  being  measured  are  lower.  This  in  and 
of  itself  creates  difficulties.  Secondly,  extensive  regions  of  laminar 
flow  may  exist  which  are  more  prone  to  separation  or  sensitivity  to 
surface  conditions  than  turbulent  boundary  layers.  Thirdly,  the 
transition  region  seems  to  be  much  longer  than  it  is  at  higher  Reynolds 
number,  and  there  is  inherent  uncertainty  introduced  by  these 
transition  effects.  As  Mueller  has  pointed  out’0,  the  free  stream 
disturbances  can  have  an  unduly  large  influence  on  the  accuracy  and 
repeatability  of  data  measured  at  low  Reynolds  numbers.  Finally,  there 
is  practically  no  dynamic  stability  data  at  all.  Because  of  the 
problems  of  taking  steady  state  data,  and  because  of  the  unknown 
effects  of  unsteady  flow  on  these  stability  derivatives,  estimated 
values  must  be  treated  with  care.  These  effects  include  the 
possibility  that  one  cannot  uncouple  the  lateral  motion  from  the 
longitudinal  motion  in  performing  the  stability  analysis,  particularly 


when  non-linear  aerodynamic  characteristics  need  to  be  used. 
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Conclusions 

Not  unexpectedly  the  usefulness  of  the  linearized  stability 
analysis  is  limited  when  applied  to  flight  at  high  lift  coefficients  at 
low  Reynolds  numbers.  The  concept  of  local  linearizations  may  well  be 
useful  as  long  as  6a  is  limited,  that  is,  it  be  the  matter  of 
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Where  these  conditions  are  violated,  then  the  linear  theory 
overestimates  both  the  frequency  of  the  motion  (i.e.,  the  period  is 
too  short)  and  the  damping  (i.e.,  6  is  too  large)  when  compared  to 
solutions  of  the  equations  of  motions.  Further  in  the  non-linear  lift 
regime  lateral  motions  will  feed  into  the  longitudinal  motion.  The 
converse  is  not  as  evident.  The  hysteresis  in  the  lift  curve  can 
couple  into  the  motion  unfavorably.  The  local  increase  in  stiffness 
coupled  with  a  drastic  reduction  in  damping  could  lead  to  divergent 
motions  in  pitch. 

The  use  of  a  generalized  phase  plane  is  indicative  of  an 
instability.  However  this  technique  is  no  more  efficient  than  any 
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trajectory  integration  scheme.  It's  virtue  is  its  compactness  of  the 
results  it  presents. 

The  clearest  need  from  this  study  is  that  of  additional  reliable 
data  in  this  Reynolds  number  range.  It  is  needed  in  several  forms: 

(1)  Static  Aerodynamic  Characteristics 

(2)  Dynamic  stability  derivatives  measured  with  small  amplitude  at  a 
variety  of  set  points. 

(3)  Some  trajectory  simulation  data  using  a  dynamic,  perhaps  a 
magnetic  model  suspension  system,  to  probe  for  unexpected,  non-linear 
aerodynamic  coupling. 

The  dynamic  response  of  the  LODED  vehicle  was  modeled  three 
different  ways.  All  methods  that  indicated  the  lateral  resonse  were 
unstable.  Liapunov  stability  analysis  proved  to  be  a  useful  way  of 
determining  the  stability  (or  instability)  of  the  system  for  both  linear 
and  nonlinear  models.  Strip  theory  with  a  wind  tunnel  correction  factor 
was  successfully  used  to  calculate  external  forces  and  moments  for  the 
full  nonlinear  control  fixed  equations.  Further  investigation  is  needed 
to  resolve  the  longitudinal/lateral  coupling,  to  fully  determine  the 
effects  of  a  hysteresis  loop  and  also  to  completely  analyze  the  results. 
The  aerodynamic  data  should  be  further  correlated  to  provide 
compatability  among  the  three  methods  of  investigations. 
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APPENDIX  A  (Ref.  7) 


Vehicle:  Loded  GW- 140  LBS,  CG-19.8  inches  aft  of  fwd  pivot  -  CGI 

Lower  winglets  with  wing  pivots  at  .32  chord  (calc.  &  meas.  der) 


Pertinent  Airplane  Characteristics 


Density  (slugs/f t**3 ) 

= 

0.00237 

Velocity  (ft/sec) 

71 .000 

Mass  (slugs) 

= 

4.350 

IYY  (slug-f t**2) 

= 

13.89 

Thrust  (pounds) 

= 

0.000 

ZJ  (ft) 

= 

0.000 

GCOS  (gamma)  (ft/sec/sec) 

= 

32.200 

GSIN  (gamma)  (ft/sec/sec) 

as 

0.000 

COS  (XZ) 

= 

0.999 

SIN  (XZ) 

=3 

0.052 

Wing  Area  (ft**2) 

10.000 

Horiz.  Tail  Area  (ft**2) 

= 

10.000 

Wing  Span  (ft) 

= 

10.000 

Horiz.  Tail  Span  (ft) 

5* 

10.000 

Wing  Chord  (ft) 

= 

1  .000 

Horz.  Tail  Chord  (ft) 

=S 

1  .000 

Wing  Aspect  Ratio 

= 

10.000 

Horz.  Tail  Aspect  Ratio 

= 

10.000 

Wing  Taper  Ratio 

= 

1.000 

Horz.  Tail  Taper  Ratio 

= 

1.000 

Wing  Alpha  (degrees) 

= 

21.800 

Tail  Alpha  (degrees) 

= 

6.667 

IWing  (degrees) 

= 

16.000 

ITail  (degrees) 

- 

6.000 

Downwash  Angle  (degrees) 

= 

5.133 

Downwash/Alpha 

= 

0.127 

Elevator  Angle  (degrees) 

-5.306 

Elevator  Area  (ft**2) 

= 

4.800 

Tail  Efficiency 

= 

1.000 

Elevator  Chord  (ft) 

= 

1.000 

2-D  Wing  CLA 

= 

0.040 

2-D  Tail  CLA 

= 

0.113 

CDO 

= 

0.050 

2-D  Wing  CDA 

= 

0.191 

2-D  Wing  CL 

z 

1.775 

Wing  CMAC 

= 

-0.076 

Distances 


Length  of  Fuselage  (ft) 

Width  of  Fuselage  (ft) 

C.G.  to  Tail  Quarter  Chord  (ft) 
Wing  to  Tail  Quarter  -  Chord  (ft) 
C.G.  to  Wing  A.C.  (chordwise)  (ft) 
C.G.  to  Wing  A.C.  (vertical)  (ft) 
Nose  to  Wing  Quarter  -  Chord  (ft) 
C.G.  to  Wing  Quarter  -  Chord  (ft) 
C.G.  to  Thrust  Axis  (ft) 


7.000 
1.000 
2.470 
4.167 
1.330 
-0.410 
1.043 
-1 .690 
0.0000 


Longitudinal  Stability  Derivatives 


CL 

= 

2.3000 

CLA 

= 

4.0110 

CLAD 

= 

5.4581 

CLQ  = 

18.3339 

CD 

JS 

0.2320 

CDA 

1.1170 

CDAD 

& 

0.0000 

CDQ  = 

0.0000 

CM 

* 

0.0000 

CMA 

= 

-9.3910 

CMAD 

=S 

-13.4815 

CMQ  = 

-74.8187 

CLDE 

= 

2.7690 

CLU 

= 

0.0000 

CT 

= 

0.0000 

ODDE 

0.7000 

CDU 

0.0000 

CTU 

3 

0.0000 

CMDE 

= 

4.6100 

CMU 

= 

0.0000 

CTRPM 

= 

0.0000 

L 
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4.79006 

0.76236 

52.99 

34.48468 

-0.90265 

-51.71737 

5.74807 

0.91483 

60.52 

35.63810 

-1.37764 

-78.93292 

6.38675 

1  .01648 

58.59 

35.35603 

-1 .76333 

-101.03116 

7.02542 

1.11813 

50.47 

34.06121 

-2.11651 

-121 .26698 

7.98344 

1.27060 

36.55 

31.25798 

-2.49515 

-142.96154 

10.00000 

1  .59155 

19.49 

25.79523 

-2.89862 

-166.07883 

100.00000 

15.91549 

0.40 

-7.95587 

-4.34941 

-249.20278 

1000.00000 

152.15494 

0.04 

-28.45762 

-4.67488 

-267.85109 

Vehicle:  Loded 

GW- 140  lbs, 

CG-19.8 

inches  aft  of 

fwd  pivot  - 

CGI 

Lower 

winglets  with  wing 

pivots  at  .32 

chord  (calc. 

&  meas.  der) 

Pertinent  Airplane  Characteristics 

Rho  -  0.00237  Wing  Area  =  10.00  Mass  =  4.35  GCOS  (gamma)  =  32.20 


U 

=  71.00 

Chord 

=  1. 

000  Span  =  10. 

000  GSIN  (gamma 

IXX 

=  10.1 

IXZ  = 

0.54 

IZZ 

22.6 

CLW 

= 

1.4792 

SA 

=  0.000 

DIH  = 

0.000 

ZC 

= 

0.450 

FUSVOL 

= 

7.417 

H 

=  1.182 

SV 

1.275 

BV 

1.670 

R1 

= 

0.569 

TR 

-  1.000 

zv 

0.055 

ETAV 

= 

1.000 

SBS 

= 

7.338 

LF 

=  7.000 

LTV  « 

2.544 

XM 

= 

2.775 

HI 

= 

1  .000 

H2 

=  1.000 

WFUS  = 

1.000 

SAH 

=3 

0.000 

CLA2DW 

- 

0.040 

BH 

=  10.000 

SH 

10.000 

TRH 

= 

1.000 

CLA2DH 

= 

0.113 

BA 

-  4.500 

CA 

1.000 

SR 

a 

0.000 

ALPHAR 

a 

5.800 

CDO 

=  0.050 

YI 

0.500 

HNOSE 

= 

1  .000 

WNOSE 

= 

1.000 

HFCY 

=  1.0000 

WFCY  = 

1  .0000 

LFCY 

a 

4.7467 

LMH  = 

4. 

7470 

HBCY 

-  1.1380 

WBCY  = 

1  .0000 

LBCY 

= 

6.0000 

Lateral  Stability  Derivatives 


CY3 

s 

-0.9740 

CLB  = 

-0.0888 

CNB  = 

0.0962 

CYP 

a 

0.1150 

CLP  = 

-1.1400 

CNP  = 

-0.1233 

CYR 

= 

0.4168 

cLr  = 

0.4040 

CNR  = 

-0.0513 

CYDA 

= 

-0.0810 

CLDA  = 

0.9773 

CNDA  = 

-0.0138 

CYDR 

= 

0.0000 

CLDR  = 

0.0000 

CNDR  = 

0.0000 

Response  to  Aileron  Deflection 

CYIN  =  -0.081000  CLIN  =  0.977300  CNIN  =  -0.013800  K  =  2  ACC  =  0.001000 

Dimensional  Stability  Derivatives 

YV  -  -0.18839  LB  =  -5.25247  NB  =  2.54455 
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APPENDIX  B 

Modeling  of  Fuselage  for  Stability  Derivatives 

The  contribution  of  the  fuselage  to  the  external  forces  and 
moments  is  required  in  the  strip  theory  method.  The  aerodynamics  of 
the  fuselage  are  developed  from  potential  flow  and  revised  for  the 
effects  of  viscosity. 

The  pressure  distribution  on  an  inclined  fuselage  by  the  method  of 
dipole  distribution  yields: 

C  ( x  i  v )  =  -2  —  (a  (x)  R2(x)) 

P  R( x )  dx 

Integration  over  v  results  in  the  lift  of  length  dx  for  the  fuselage: 

—  =  2irq™d(a(x)  R2(x))/dx 
dx 

For  non-circular  cross-sections,  the  semi-width,  bf(x)/2,  is 
substituted  for  R(x).  Thus,  after  integration 

2  If 

Lf  =  2irqoo(a(x)  bf(x)/4)|Q 
2 

L  “  (tt/2)  q®bf  a^, 

Xf  f 

where  subscript  If  indicates  evaluation  at  tail  of  fuselage. 

The  pitching  moment  is  obtained: 
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M 


F 


dL 

F 

-  -  x  dx 

dx 


(ir/2 ) 


c 


2 

q~a { x )  bf ( x ) 


dx 


For  a  first  approximation,  assume  a(x)  is  constant  along  the  fuselage 
and  equal  to  at  the  c.g. 

The  lift  and  pitching  moment  for  the  LODED  fuselage  are  obtained 
with  the  appropriate  width  distribution,  bf(x).  The  drawings  of  the 
fuselage  indicate  it  is  not  a  completely  closed  body  at  the  tail 
(exclusive  of  the  propeller)  with  Bf  =6  inches.  According  to 
Hafer9,  viscous  effects  can  be  accounted  for  by  using  the  geometric 
width  plus  the  boundary-layer  displacement  thickness,  6(x),  as  the 
actual  width  in  the  equations. 

From  these  considerations  the  fuselage  can  be  approximated: 

Thus, 

lF  =  (ir/2 )  q<x>a  ( 1  )2 
and  integration  of  Eq.  4: 

M  =  (ir/2)  q  a  [-  -1  +  4rx2J  |r  +  x  |  f 

F  cd  3  o  o 

substitution  of  proper  values: 

Mf  -  (6.58333)  (ir/2)  q^a 

At  zero  angle  of  attack,  the  drag  on  the  fuselage  is  nearly  equal 
to  the  drag  on  a  flat  plate  of  equal  surface  area  and  length.  From 


i 
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Ref.  (2) 


D 

F 


flat 

plate 


(1  + 


(.5 >«_) 

F 


(a-0°) 


with  6F 


fmax 


1  f t./7  ft. 


flat 

plate 


CfSF^ 


With  Reynolds  number  for  the  fuselage  =  285,000xlf/c  »  1,995,000,  the 
friction  factor  =  .0064. 


SF  -  2irr2  +  (If  -  r)(bfmax) 
S  =  27.5708  ft. 


According  to  Wieselberger2 ,  the  induced  drag  on  an  inclined  fuselage 
is  estimated: 


D 


2 


pV27T 


B 


<B2-b2  )‘ 

fmax 


B  -  semi  span  “  5  ft. 


Thus  the  total  drag  on  a  fuselage  at  an  angle  of  attack: 

S2 

2 


°F"  cfV1+(*5)Vq~  F 


B 


2  2  2' 

pV  w  (B  -b“  ) 

fmax 


Substitution  of  appropriate  constants  leads  to: 

.  2 


D  =■  ( . 1 888 )q  +  (3.9788736) 


oV 
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The  side  force,  Y ,  and  yawing  moment,  N,  are  dependent  upon  the 
sideslip  angle,  3,  and  are  modeled  analogously  to  the  lift  and  pitching 
moment,  respectively: 

yF  **  -(ir/2)q<o6 
nF  =»  -(6.58333)  (rr/2 )  q^a 
Equations  5  to  9  are  implemented  in  the  program. 
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TABLE  I 

RUNS  FOR  LINEARIZED  ANALYSIS 


6u,  6v,  Time  of 


RUN 

#  6a 

fiw,  needed 

Disturb 

1  o 

n 

u 

a> 

0 

II 

0° 

1 

+i° 

fiw  =  1 .286 

.1 

sec 

2 

+i° 

fiw  =  1 .286 

.5 

sec 

3 

+2° 

fiw  =  2.582 

.1 

sec 

4 

+2° 

fiw  =  2.582 

.5 

sec 

5 

+4° 

fiw  «  5.217 

.  1 

sec 

6 

+4° 

fiw  =  5.217 

.5 

sec 

II 

U 

Q  O  “ 

0° 

7 

+  i° 

fiu 

7.454 

.1 

sec 

8 

-r 

fiw 

= 

-1 .249 

.5 

sec 

9 

+2° 

fiw 

28 

2.521 

.1 

sec 

10 

+2° 

fiw 

= 

2.521 

.5 

sec 

1  1 

-4° 

fiw 

S 

-4.97 

.1 

sec 

12 

+4° 

fiu 

= 

-22.76 

.5 

sec 

Lateral 

CL  -  1 -0  0O  =  0° 

13  68  -  +1°  6v  =  1  .201  .5  sec 

CL  =  .6  80  =  0° 

14  58  =  +4°  6 v  =  4.91  .1  sec 


AMATRIX=  fit  =■  .01  sec 

fi^u,  fir.,  _fiw  Variables  Plotted 

fit  fit  fit  Short  Period  Long  Period 

aQ  =  14.285° 


128.556 

a , 

8 

u, 

q 

128.556 

a , 

8 

u, 

8 

258.244 

a  * 

8 

u, 

8 

258.244 

a, 

8 

u, 

q 

521.71 

a , 

q 

u, 

8 

521  .71 

a , 

q 

u, 

q 

aQ  =  8.571° 


-745.44 

a  , 

8 

u, 

8 

-124.955 

a  , 

q 

u, 

q 

252.111 

a  , 

e 

u. 

e 

252.111 

a, 

q 

u, 

8 

-496.804 

a  , 

q 

u, 

q 

2275.704 

a  , 

q 

u, 

q 

B0  =  0° 

120.1  6 


490.935  8,  p  8 , 
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TABLE  II 

RUNS  FOR  PITCH-PLANE  ANALYSIS 


Duration  of 

RUN  #  Disturbance  ^Disturbance  =  68/5t 


6t  =  .005  sec 


CL  -  -6 


®o 


=  0° 


1  . 1  sec 

2  . 1  sec 

5  .5  sec 

6  .5  sec 


1  .0 
5.0 
1  .0 
5.0 


3 

.1 

sec 

1  .0 

4 

.  1 

sec 

5.0 

7 

.5 

sec 

1  .0 

8 

.5 

sec 

5.0 

9 

.5 

sec 

-5.0 

10 

.1 

sec 

-5.0 

37 


TABLE  III 
STRIP  THEORY  RUNS 


RUN  # 

cL 

At 

6  ( ctcal^ 

6  (  B  ca 1 ^ 

1 

2.33 

.1 

i° 

0° 

2 

2.33 

.  1 

2° 

0° 

3 

2.33 

.1 

4° 

0° 

4 

2.33 

.5 

1° 

0° 

5 

2.33 

.5 

2° 

0° 

6 

2.33 

.5 

4° 

0° 

7 

2.77 

.1 

1° 

0° 

8 

2.77 

.1 

2° 

0° 

9 

2.77 

.1 

4° 

0° 

10 

2.77 

.5 

r 

0° 

1  1 

2.77 

.5 

2° 

0° 

12 

2.77 

.5 

4° 

0° 

13 

2.33 

.1 

0° 

2° 

14  Lissaman 

2.33 

.1 

-1° 

0° 

FROM  NOSE 


DETAILS  OF  LODED  FUSELAGE  GEOMETRY 


[G  DURATION  EXPENDABLE  DECOY 


2.8 


f-IG.  3C  COMPARISON  OF  LIFT  CURVE  FOR  STABILITY  DERIVATIVE  AND  MEASURED  LIFT  CURVE 
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FIG.  4  COORDINATE  SYSTEM 


FIG.  5  PHASE  PLANE  -  A  STABLE  SPIRAL 


SHORT  PERIOD  MOTION  (LINEARISED) 


Time  (Sec) 


FIG.  7b  PITCH  PLANE  SOLUTION  FOR  LONG  TIME  OR  PHUGOID  MOTION  (LINEARIZED) 


FIG.  8A  NON-LINEAR  MOTION  C  =  2.33,  6=1®  ANGLE  OF  ATTACK  VARIATION 


-0.8 
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.33,  5=1°  ANGULAR  VELOCITY  VARIATION 


2.0 


.33,  <5  =  1°  LONG  TIME  VELOCITY  VARIATION 


Ve  locit 


2.0 


FIG.  11  LIAPUNOV  STABILITY  C  =2.33 


0.5 
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FIG.  12  LIAPUNOV  STABILITY 


FIG.  14  LIAPUNOV  STABILITY  C  =2.77 


Velocity 
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f 


0.25 


0.20 


0.15 


0.10 


0.05 


O 


-0.05 


-0.10 
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FIG.  15  LIAPUNOV  STABILITY  =  2.33,  (5  =  0°.  YAW  =  2° 


